Smart-seq2 Quality Control

Gametocyte Development in Plasmodium berghei

---

1. Read in the data

Load the required packages

[1] "Seurat is loaded correctly"
[1] "viridis is loaded correctly"
[1] "ggpubr is loaded correctly"
[1] "cowplot is loaded correctly"
[1] "ggplot2bdc is loaded correctly"
[1] "patchwork is loaded correctly"

Import GTF file

This will be helpful later on. This contains annotations for each gene:

##Import gtf file:
gtf <- read.table("/Users/Andy/GCSKO/GCSKO_analysis_git/data/Reference/Pberghei.gtf", sep="\t", header = FALSE)

## inspect
head(gtf)

*N.B.This was the .gtf file used for feature counting the data

Import Counts and Pheno files

## read in counts
counts <- read.delim("/Users/Andy/GCSKO/GCSKO_analysis_git/data/Smartseq2/counts_2020.csv", sep = ",", stringsAsFactors = FALSE, row.names = 1)

## read in pheno
pheno <- read.delim("/Users/Andy/GCSKO/GCSKO_analysis_git/data/Smartseq2/pheno_2020.csv", sep = ",", stringsAsFactors = FALSE)

## check dimensions of both
expected_number_of_cells <- ( 384 * 12)
paste("The expected number of cells is", expected_number_of_cells)

[1] "The expected number of cells is 4608"

paste("number of genes in counts", dim(counts)[1], "number of cells in counts", dim(counts)[2])

[1] "number of genes in counts 5250 number of cells in counts 4608"

paste("number of columns in counts", dim(pheno)[2], "number of cells in pheno", dim(pheno)[1])

[1] "number of columns in counts 98 number of cells in pheno 4608"

## make rownames the sample_id in pheno
rownames(pheno) <- pheno$sample_id

## check that the names are the same
paste("Are the names the same?")

[1] "Are the names the same?"

table((colnames(counts) %in% rownames(pheno)))

TRUE 
4608 

Remove Anopheles cells

These cells were included in one of the sequencing runs and are irrelevant.

## get rid of the cells by finding all the cells that have As in their species col
names_of_mosquito_cells <- rownames(pheno[(pheno$species == "As"), ])

## which gives
table(colnames(counts) %in% names_of_mosquito_cells)

FALSE  TRUE 
 4512    96 

## subset 
counts <- counts[,-which(colnames(counts) %in% names_of_mosquito_cells)]
pheno <- pheno[-(which(rownames(pheno) %in% names_of_mosquito_cells)), ]

## check dimensions
paste("number of genes in counts", dim(counts)[1], "number of cells in counts", dim(counts)[2])

[1] "number of genes in counts 5250 number of cells in counts 4512"

paste("number of columns in counts", dim(pheno)[2], "number of cells in pheno", dim(pheno)[1])

[1] "number of columns in counts 98 number of cells in pheno 4512"

Remove control cells:

## find all cells that are labelled as such
names_of_control_cells <- rownames(pheno[(pheno$is_control == "TRUE"), ])

## which gives
table(colnames(counts) %in% names_of_control_cells)

FALSE  TRUE 
 4386   126 

## subset 
counts <- counts[,-which(colnames(counts) %in% names_of_control_cells)]
pheno <- pheno[-(which(rownames(pheno) %in% names_of_control_cells)), ]

## check dimensions
paste("number of genes in counts", dim(counts)[1], "number of cells in counts", dim(counts)[2])

[1] "number of genes in counts 5250 number of cells in counts 4386"

paste("number of columns in counts", dim(pheno)[2], "number of cells in pheno", dim(pheno)[1])

[1] "number of columns in counts 98 number of cells in pheno 4386"

Remove 10 and 27

10 and 27 failed genotyping and so will be removed

## find all cells that are labelled as such
names_of_genotype_fail_cells <- rownames(pheno[(pheno$sub_identity_updated %in% c("GCSKO-10", "WT-10", "GCSKO-27")), ])

## which gives
table(colnames(counts) %in% names_of_genotype_fail_cells)

FALSE  TRUE 
 3732   654 

## subset 
counts <- counts[,-which(colnames(counts) %in% names_of_genotype_fail_cells)]
pheno <- pheno[-(which(rownames(pheno) %in% names_of_genotype_fail_cells)), ]

## check dimensions
paste("number of genes in counts", dim(counts)[1], "number of cells in counts", dim(counts)[2])

[1] "number of genes in counts 5250 number of cells in counts 3732"

paste("number of columns in counts", dim(pheno)[2], "number of cells in pheno", dim(pheno)[1])

[1] "number of columns in counts 98 number of cells in pheno 3732"

Filter non-genes out of counts

Filter out the QC columns contained in counts that do not correspond to genes

## since the datatable currently contains:
paste("The current counts table contains these non-gene rows:")

[1] "The current counts table contains these non-gene rows:"

rownames(counts[-(grep("PBANKA*", rownames(counts))), ])

[1] "__alignment_not_unique" "__ambiguous"            "__no_feature"          
[4] "__not_aligned"          "__too_low_aQual"       

## subset
# make a copy of counts just in case that info is useful in future analysis
counts_genes <- counts
# remove non-gene rows
counts <- counts[(grep("PBANKA*", rownames(counts))), ]
#check
paste("number of genes in counts", dim(counts)[1], "number of cells in counts", dim(counts)[2])

[1] "number of genes in counts 5245 number of cells in counts 3732"

This 5245 includes all types of genes:

## check size
paste("The dimensions of the GTF file are:")

[1] "The dimensions of the GTF file are:"

dim(gtf[-which(gtf$V3 == "CDS"),])

[1] 5245    9

## look at possible types of genes
paste("The types of genes in the dataframe are:")

[1] "The types of genes in the dataframe are:"

names(table(gtf$V3))

[1] "CDS"                    "mRNA"                   "ncRNA"                 
[4] "pseudogenic_transcript" "rRNA"                   "snoRNA"                
[7] "snRNA"                  "tRNA"                  

---

2. Filtering

set-up data

make seurat object

## make seurat object
ss2_mutants <- CreateSeuratObject(counts = counts, project = "GCSKO", min.cells = 0, min.features = 1, meta.data = pheno)

Feature names cannot have underscores ('_'), replacing with dashes ('-')

## n.b. cells must have at least 1 gene per cell because otherwise this causes problems downstream

## add experiment column to meta data:
ss2_mutants@meta.data$experiment <- "ss2_mutants"

## inspect object
ss2_mutants

An object of class Seurat 
5245 features across 3705 samples within 1 assay 
Active assay: RNA (5245 features, 0 variable features)

## how many genes are not detected at all?
paste(length(which(rowSums(as.matrix(ss2_mutants@assays$RNA@counts)) == 0)), "genes are not detected in any cell")

[1] "63 genes are not detected in any cell"

[1] "27 cells have zero counts"

Filtered Cells Plate Heatmap

set up data

## make dataframe
df_platemap <- data.frame(well_name = pheno$well_name_R, cell_name = rownames(pheno), row.names = rownames(pheno))

## get cells that are filtered out
cells_kept <- which(rownames(df_platemap) %in% rownames(ss2_mutants@meta.data))

## make extra column in plotting df
df_platemap$colour <- "filtered_out"
df_platemap$colour[cells_kept] <- "passed_QC"

## inspect
#head(df_platemap)

##Make a table for each

##make a new dataframe for plotting
table_platemap <- as.data.frame(table(df_platemap$well_name, df_platemap$colour))

## modify names of dataframe
names(table_platemap) <- c("well", "filtered", "frequency")

#make separate dfs
df_1 <- table_platemap[table_platemap$filtered == "filtered_out",]
df_2 <- table_platemap[table_platemap$filtered == "passed_QC",]

# merge them back together
table_platemap <- merge(df_1, df_2, by = "well", all = TRUE)

## remove the irrelevant rows and rename cols and rename rows
table_platemap <- table_platemap[ ,c(1,3,5)]
names(table_platemap) <- c("well", "filtered_out", "passed_QC")
table_platemap$well_name <- rownames(table_platemap)

## add A1?
#table_platemap[96,] <- c("A1", "0", "0")

## add cols for plotting
table_platemap$Row <- as.numeric(match(toupper(substr(table_platemap$well, 1, 1)), LETTERS))
table_platemap$Column <- as.numeric(substr(table_platemap$well, 2, 5))

## calculate a percentage
table_platemap$filtered_pc <- ((table_platemap$filtered_out)/(table_platemap$filtered_out + table_platemap$passed_QC))*100

## due to an error in the scale_y_reverse, you need to make a column of 'reversed' row values to plot
for(i in seq(1,8)){
  table_platemap$Row_rev[table_platemap$Row == i] <- seq(8,1)[i]
}

## find wells where it's zero
zero_wells <- table_platemap$filtered_pc == 0
table_platemap$filtered_pc[zero_wells] <- NA

Plot

## plot
sample_map_no_reads <- ggplot(data=table_platemap, aes(x=Column, y=Row_rev)) +
  #set up the platemap layout
    geom_point(data=expand.grid(seq(1,12), seq(1,8)), aes(x=Var1, y=Var2), color="grey90", fill="white", shape=21, size=8) +
  #Change the shape and colour of points for a variable
    geom_point(aes(colour = filtered_pc), size = 7) +
  #change the colours
    scale_colour_viridis_c(guide = "colourbar", na.value="white") +
  #fix the ratio of coordinates
    coord_fixed(ratio=(13/12)/(9/8), xlim=c(0.5, 12.5), ylim=c(0.6, 8.4)) +
  # make into a plate plot
    theme_bdc_microtiter() +
  #add labels for the y axis
    scale_y_continuous(breaks = seq(1, 8), labels = LETTERS[8:1]) +
  #add labels for the x axis
    scale_x_continuous(position = "top", breaks=seq(1, 12)) +
  #Add a title and change label of fill
    labs(title="The position of cells with zero counts" , size = 6, colour = "percentage of cells in this well with zero counts") +
  #Improve aes of guides
    guides(colour = guide_colorbar(barwidth = 4, barheight = 1, title="percentage of cells in this well with zero counts")) +
# display legend on bottom and make text bigger
    theme(legend.position="bottom", legend.text=element_text(size=10), legend.title = element_text(size = 10))

## print
sample_map_no_reads

## save
ggsave("/Users/Andy/GCSKO/GCSKO_analysis_git/images_to_export/QC_SS2_platemap_no_count.png", plot = sample_map_no_reads, device = "png", path = NULL, scale = 1, width = 12, height = 12, units = "cm", dpi = 300, limitsize = TRUE)

Filter Mitochondrial

Mitochondrial cell counts

## old method
## This old method just picked out rRNA genes as a rough estimate, the method below picks out all genes expressed from the mitchondrial genome

## extract mitochondrial genes 
#mito_genes <- gtf[which(gtf$V3 == "rRNA"),]$V9
#mito_genes <- gsub(";.*","", gsub("gene_id ", "", mito_genes))
#paste("These are the mitochondrial genes")
#head(mito_genes)

## make a percentage mitocondrial for each cell (this will work as long as you filter cells out with zero counts)
#ss2_mutants <- PercentageFeatureSet(ss2_mutants, features = which(rownames(counts) %in% mito_genes), col.name = "percent.mt")

## extract mito genes
mito_genes <- ss2_mutants@assays$RNA@counts@Dimnames[[1]][grep("^PBANKA-MIT",ss2_mutants@assays$RNA@counts@Dimnames[[1]])]

## plot the genes individually
VlnPlot(object = ss2_mutants, features = mito_genes, pt.size = 0.01, group.by = "experiment")

All cells have the same value of PBANKA-MIT00100:rRNA.All cells have the same value of PBANKA-MIT00200:rRNA.All cells have the same value of PBANKA-MIT00900:rRNA.All cells have the same value of PBANKA-MIT02000:rRNA.All cells have the same value of PBANKA-MIT02600:rRNA.All cells have the same value of PBANKA-MIT03000:rRNA.All cells have the same value of PBANKA-MIT03300:rRNA.All cells have the same value of PBANKA-MIT03400:rRNA.

## make a percentage mitocondrial for each cell (this will work as long as you filter cells out with zero counts)
ss2_mutants <- PercentageFeatureSet(ss2_mutants, pattern = "^PBANKA-MIT", col.name = "percent.mt")

plot percentage mitochondrial

## plot for percentage of mitochondrial reads
v1 <- VlnPlot(object = ss2_mutants, features = "percent.mt", group.by = "sub_identity_updated", pt.size = 0.01) + 
  ## add a line where we will filter 
  geom_hline(yintercept=20) +
  ## remove legend
  theme(legend.position = "none") +
  ## change labels
  labs(x = "Genotype", y = "% Mitochondrial Reads") +
  ## remove plot title
  theme(plot.title = element_blank())

## plot for percentage of mitochondrial reads
v2 <- VlnPlot(object = ss2_mutants, features = "percent.mt", group.by = "experiment", pt.size = 0.01) +
  geom_hline(yintercept=20) +
  ## remove legend
  theme(legend.position = "none") +
  ## fremove axis labels
  labs(x="", y = "") +
  ## remove axis elements
  theme(plot.title = element_blank(), axis.text.y = element_blank(),
        axis.ticks.y = element_blank()) +
  ## change label on bottom of plot
  scale_x_discrete(labels = "All cells")

## plot together
QC_mito_violin <- v1 + v2 + plot_layout(ncol = 2, nrow = 1, widths = c(4, 1), heights = c(2, 2))

## print
QC_mito_violin

## save
ggsave("/Users/Andy/GCSKO/GCSKO_analysis_git/images_to_export/QC_SS2_mito_violin.png", plot = QC_mito_violin, device = "png", path = NULL, scale = 1, width = 15, height = 10, units = "cm", dpi = 300, limitsize = TRUE)

make a dataframe for plotting

df <- data.frame(nCount = log10(ss2_mutants@meta.data$nCount_RNA), nGenes = ss2_mutants@meta.data$nFeature_RNA, identity = ss2_mutants@meta.data$identity_updated, identity_gene = ss2_mutants@meta.data$identity_gene_updated, percent_mt = ss2_mutants@meta.data$percent.mt)

Where do mitchondrial poor cells lie?

## make extra col for filtered values:
filtered_out_cells <- which(df$percent_mt > 20) 

## plot
QC_mito_graph <- ggplot(df, aes(x = nCount, y = nGenes)) + 
  geom_point(alpha = 0.4) +
  scale_size(range = c(0.1,4)) +
  #geom_rug() + 
  scale_y_continuous(name = "Number of Detected Genes") + 
  scale_x_continuous(name = "log10(Number of Total Counts)") +
  scale_color_viridis(option = "D") +
  theme_pubr() +
  theme(legend.position = "bottom") +
  geom_vline(xintercept=3) +
  geom_hline(yintercept=220) +
  geom_hline(yintercept=3300) +
  #annotate selected points
  annotate("point", df$nCount[filtered_out_cells], df$nGenes[filtered_out_cells], size = 2, colour = "orange")

## add this to gemo_point if needed: aes(colour = percent_mt, size = percent_mt)

## print
QC_mito_graph

## save
ggsave("/Users/Andy/GCSKO/GCSKO_analysis_git/images_to_export/QC_SS2_mito_graph.png", plot = QC_mito_graph, device = "png", path = NULL, scale = 1, width = 10, height = 10, units = "cm", dpi = 300, limitsize = TRUE)

[1] "The nGene and nCounts filters catch most of the high % mito cells ( 115 ) but introducing a filter would capture a further 68"

## plot 1
plot1 <- FeatureScatter(ss2_mutants, feature1 = "nCount_RNA", feature2 = "percent.mt", pt.size = 0.01, group.by = "identity_updated")

## plot 2
plot2 <- FeatureScatter(ss2_mutants, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", pt.size = 0.01, group.by = "identity_updated")

## plot together
plot1 + plot2

Filter poor quality cells

in the MCA paper: “Cells with fewer than 1000 genes per cell and 2500 reads per cell were removed from the liver-stage parasites, trophozoites, male and female gametocytes, ookinetes, ookinetes/oocysts, and oocyst stages. Cells with fewer than 500 genes per cell and 2500 reads per cell were removed from schizonts and injected sporozoites. Cells with fewer than 40 genes per cell and 1000 reads per cell were removed from merozoites, rings, and gland sporozoites (fig. S2 and table S1). Additionally, we removed genes from further analysis that were detected in fewer than two cells across the entire dataset. The final dataset contained 1787 high-quality single cells from 1982 sequenced cells and 5156 genes out of 5245 genes with annotated transcripts.”

so use: 1000 reads per cell & 40 genes per cell as absolute minimums

## plot main dotplot
plot1 <- ggplot(df, aes(x = nCount, y = nGenes, color = identity)) + 
  geom_point(aes(color = identity), size = 0.1) +
  geom_rug() + 
  scale_y_continuous(name = "Number of Detected Genes") + 
  scale_x_continuous(name = "log10(Number of Total Counts)") + 
  theme_pubr() +
  theme(legend.position = "bottom") +
  geom_vline(xintercept=3) +
  geom_hline(yintercept=220) +
  geom_hline(yintercept=3300) +
  ## improve aesthetics of guide
  guides(colour = guide_legend(override.aes = list(size=4)), title="Genotype")

## plot density plot x
dens1 <- ggplot(df, aes(x = nCount, fill = identity)) + 
  geom_density(alpha = 0.2) + 
  theme_void() + 
  theme(legend.position = "none")

## plot density plot y
dens2 <- ggplot(df, aes(x = nGenes, fill = identity)) + 
  geom_density(alpha = 0.2) + 
  theme_void() + 
  theme(legend.position = "none") + 
  coord_flip()

## plot together
QC_composite_plot <- dens1 + plot_spacer() + plot1 + dens2 + plot_layout(ncol = 2, nrow = 2, widths = c(4, 1), heights = c(1, 4))

## print
QC_composite_plot

## save
ggsave("/Users/Andy/GCSKO/GCSKO_analysis_git/images_to_export/QC_SS2_composite_plot.png", plot = QC_composite_plot, device = "png", path = NULL, scale = 1, width = 20, height = 20, units = "cm", dpi = 300, limitsize = TRUE)

---

3. Analysis

Final filter

filter out cells

## subset
ss2_mutants_final <- subset(ss2_mutants, subset = nFeature_RNA > 220 & nFeature_RNA < 3300 & percent.mt < 20 & nCount_RNA > 1000)

## inspect resultant object
ss2_mutants_final

An object of class Seurat 
5245 features across 2970 samples within 1 assay 
Active assay: RNA (5245 features, 0 variable features)

[1] "735  cells are removed by these filters"

Platemap heatmap for failed cells

set up data

## make dataframe
df_platemap <- data.frame(well_name = pheno$well_name_R, cell_name = rownames(pheno), row.names = rownames(pheno))

## get cells that are filtered out
cells_kept <- which(rownames(df_platemap) %in% rownames(ss2_mutants_final@meta.data))

## make extra column in plotting df
df_platemap$colour <- "filtered_out"
df_platemap$colour[cells_kept] <- "passed_QC"

## inspect
head(df_platemap)

## make a new dataframe with metrics of interest
table_platemap <- as.data.frame(table(df_platemap$well_name, df_platemap$colour))
names(table_platemap) <- c("well", "filtered", "frequency")

# make separate dfs
df_1 <- table_platemap[table_platemap$filtered == "filtered_out",]
df_2 <- table_platemap[table_platemap$filtered == "passed_QC",]

# merge them back together
table_platemap <- merge(df_1, df_2, by = "well", all = TRUE)

## remove the irrelevant rows and rename cols and rename rows
table_platemap <- table_platemap[ ,c(1,3,5)]
names(table_platemap) <- c("well", "filtered_out", "passed_QC")
table_platemap$well_name <- rownames(table_platemap)

## add cols for plotting
table_platemap$Row <- as.numeric(match(toupper(substr(table_platemap$well, 1, 1)), LETTERS))
table_platemap$Column <- as.numeric(substr(table_platemap$well, 2, 5))

## inspect:
head(table_platemap)

## calculate percentage
table_platemap$filtered_pc <- ((table_platemap$filtered_out)/(table_platemap$filtered_out + table_platemap$passed_QC))*100

## due to an error in the scale_y_reverse, you need to make a column of 'reversed' row values to plot
for(i in seq(1,8)){
  table_platemap$Row_rev[table_platemap$Row == i] <- seq(8,1)[i]
}

zero_wells <- table_platemap$filtered_pc == 0
table_platemap$filtered_pc[zero_wells] <- NA

plot

## plot
platemap_failed_qc <- ggplot(data=table_platemap, aes(x=Column, y=Row_rev)) +
  #set up the platemap layout
    geom_point(data=expand.grid(seq(1,12), seq(1,8)), aes(x=Var1, y=Var2), color="grey90", fill="white", shape=21, size=8) +
  #Change the shape and colour of points for a variable
    geom_point(aes(colour = filtered_pc), size = 7) +
  #change the colours
    scale_colour_viridis_c(guide = "colourbar", na.value="white") +
  #Add a title and change label of fill
    labs(title="The position of cells that failed QC", size = 6) +
  #fix the ratio of coordinates
    coord_fixed(ratio=(13/12)/(9/8), xlim=c(0.5, 12.5), ylim=c(0.6, 8.4)) +
  # make into a plate plot
    theme_bdc_microtiter() +
  #add labels for the y axis
    scale_y_continuous(breaks = seq(1, 8), labels = LETTERS[8:1]) +
  #add labels for the x axis
    scale_x_continuous(position = "top", breaks=seq(1, 12)) +
  #Improve aes of guides
    guides(colour = guide_colorbar(barwidth = 4, barheight = 1, title="percentage of cells in well that failed QC")) +
  # display legend on bottom and make text bigger
    theme(legend.position="bottom", legend.text=element_text(size=10), legend.title = element_text(size = 9))

## print
platemap_failed_qc

## save
ggsave("/Users/Andy/GCSKO/GCSKO_analysis_git/images_to_export/QC_SS2_platemap_failed.png", plot = platemap_failed_qc, device = "png", path = NULL, scale = 1, width = 10, height = 10, units = "cm", dpi = 300, limitsize = TRUE)

bar chart of distance from edge of plate

Make a bar graph that describes the number of failed cells by distance from the edge of the plate

## get lists of cells for each square around the plate
total_0 <- table_platemap[which(table_platemap$Row %in% c(1,8) | table_platemap$Column %in% c(1,12)), ]

total_1 <- table_platemap[which(table_platemap$Row %in% c(2,7) | table_platemap$Column %in% c(2,11)), ]
total_1 <- total_1[-which(total_1$well %in% total_0$well),]

total_2 <- table_platemap[which(table_platemap$Row %in% c(3,6) | table_platemap$Column %in% c(3,10)), ]
total_2 <- total_2[-which(total_2$well %in% total_0$well | total_2$well %in% total_1$well),]

total_3 <- table_platemap[which(table_platemap$Row %in% c(4,5) | table_platemap$Column %in% c(4,9)), ]
total_3 <- total_3[-which(total_3$well %in% total_0$well | total_3$well %in% total_1$well | total_3$well %in% total_2$well),]

## create a dataframe which calculates the percentage failed for each distance from the edge
frequency_of_failed_cells <- data.frame(percentage_failed = 
                                          c(
                                            sum(total_0$filtered_out) / (sum(total_0$passed_QC) + sum(total_0$filtered_out)),
                                            sum(total_1$filtered_out) / (sum(total_1$passed_QC) + sum(total_1$filtered_out)),
                                            sum(total_2$filtered_out) / (sum(total_2$passed_QC) + sum(total_2$filtered_out)),
                                            sum(total_3$filtered_out) / (sum(total_3$passed_QC) + sum(total_3$filtered_out))
                                          ))

## neaten the dataframe
frequency_of_failed_cells$distance_from_edge <- c(0,1,2,3)

position_qc_fail <- ggplot(frequency_of_failed_cells, aes(x=distance_from_edge, y=percentage_failed)) +
  geom_bar(stat="identity", fill="black")+
  theme_classic() +
  theme(panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), text = element_text(size=12)) +
  ## change the axis titles
  labs(y= "Cells that failed QC (%)", x = "Number of wells from the edge of the plate") +
  ## prevent the title going off the edge of the plot
  coord_cartesian(clip = 'off')

## save
ggsave("/Users/Andy/GCSKO/GCSKO_analysis_git/images_to_export/QC_SS2_failed_position.png", plot = position_qc_fail, device = "png", path = NULL, scale = 1, width = 10, height = 10, units = "cm", dpi = 300, limitsize = TRUE)

## plot
position_qc_fail

Mapping Plot (Figure S6A)

prepare

## make dataframe
df_mapping_rates_hisat <- data.frame(alignment_rate = pheno$overall_hisat2_alignment_rate, cell_name = rownames(pheno), row.names = rownames(pheno))

## get cells that are filtered out
cells_filtered <- which(rownames(df_mapping_rates_hisat) %in% rownames(ss2_mutants_final@meta.data))

## make extra column in plotting df
df_mapping_rates_hisat$colour <- "Failed QC"
df_mapping_rates_hisat$colour[cells_filtered] <- "Passed QC"

## plot
mapping_rate_plot <- ggplot(df_mapping_rates_hisat, aes(x = reorder(cell_name,alignment_rate,sum), y=alignment_rate, fill = colour)) +
  geom_col() + 
  theme(legend.position = "none") +
  labs(x="Cells", y = "Overall Alignment Rate (%)", size = 15, fill = "QC Result") +
  theme_classic() +
  theme(plot.title = element_blank(), axis.text.x = element_blank(),
        axis.ticks.x = element_blank(), text = element_text(size=20)) +
  ## change the colours
   scale_fill_manual(values=c("red2", "green4"))

## save
ggsave("/Users/Andy/GCSKO/GCSKO_analysis_git/images_to_export/QC_SS2_mapping_rate.eps", plot = mapping_rate_plot, device = "eps", path = NULL, scale = 1, width = 20, height = 10, units = "cm", dpi = 300, limitsize = TRUE)

## plot
mapping_rate_plot

Recovery plots

How many cells are recovered per condition?

## extract data
df_cells_recovered <- as.data.frame(table(ss2_mutants_final@meta.data$sub_identity_updated))
df_cells_profiled <- as.data.frame(table(pheno$sub_identity_updated))

## make a df of the metrics
df_cells_recovered <- merge(df_cells_recovered, df_cells_profiled, by = "Var1")
names(df_cells_recovered) <- c("genotype", "recovered", "profiled")

## calculate perecentage
df_cells_recovered$recovered_pc <- (df_cells_recovered$recovered/df_cells_recovered$profiled)*100

## round up percentage for labels
df_cells_recovered$recovered_pc_rounded <- signif(df_cells_recovered$recovered_pc, 2)
df_cells_recovered$recovered_pc_rounded <- paste0(df_cells_recovered$recovered_pc_rounded, " %")

## make a column for number of cells annotation
df_cells_recovered$number_cells_annotation <- paste0("n = ", df_cells_recovered$recovered)

## inspect dataframe
df_cells_recovered

##plot
QC_by_genotype <- ggplot(df_cells_recovered, aes(x=genotype, y=recovered_pc ,fill=genotype)) + 
  geom_bar(stat = "identity") +
  ylim(0, 100) +
  #scale_fill_brewer(palette = "Set1") +
  geom_text(aes(label=recovered_pc_rounded, vjust=0.5, hjust = 1.5)) +
  geom_text(aes(label=number_cells_annotation, vjust=0.5, hjust = -0.05)) +
  labs(y = "Percentage of cells recovered", x = "Genotype") +
  coord_flip(clip = "off") +
  theme_classic() +
  ## remove legend + add border so it doesn't clip the annotation 
  theme(legend.position="none", plot.margin = unit(c(1,3,1,3), "lines"))

## plot
QC_by_genotype

## save
ggsave("/Users/Andy/GCSKO/GCSKO_analysis_git/images_to_export/QC_SS2_by_genotype.png", plot = QC_by_genotype, device = "png", path = NULL, scale = 1, width = 15, height = 10, units = "cm", dpi = 300, limitsize = TRUE)

Investigate 20 and 2 further

## how many plates?
plate_investigation <- pheno[which(pheno$sub_identity_updated == "GCSKO-20" | pheno$sub_identity_updated == "GCSKO-2" | pheno$sub_identity_updated == "WT-2" | pheno$sub_identity_updated == "WT-20"), ]
table(plate_investigation$plate_id_unique)

GCSKO-10_plate_1 GCSKO-10_plate_2 GCSKO-10_plate_3 GCSKO-10_plate_4 
              93               93               93               93 
 GCSKO-2_plate_1  GCSKO-2_plate_2 
              92               92 

add columns if it failed QC:

## get cells that are filtered out
cells_kept <- which(rownames(plate_investigation) %in% rownames(ss2_mutants_final@meta.data))

## make extra column in plotting df
plate_investigation$colour <- "filtered_out"
plate_investigation$colour[cells_kept] <- "passed_QC"

table(plate_investigation$colour)

filtered_out    passed_QC 
         359          197 

## add cols for plotting
plate_investigation$Row <- as.numeric(match(toupper(substr(plate_investigation$well_name_R, 1, 1)), LETTERS))
plate_investigation$Column <- as.numeric(substr(plate_investigation$well_name_R, 2, 5))

plot

## find unique plate names
plate_names <- unique(plate_investigation$plate_id_unique)

## due to an error in the scale_y_reverse, you need to make a column of 'reversed' row values to plot
for(i in seq(1,8)){
  plate_investigation$Row_rev[plate_investigation$Row == i] <- seq(8,1)[i]
}

## make an empty list
plot_list = list()

## run for loop which will generate plots for each plate
for(i in seq_along(plate_names)){
## make df
  table_platemap <- plate_investigation[plate_investigation$plate_id_unique == plate_names[i], ]
  ## plot
  p <- ggplot(data=table_platemap, aes(x=Column, y=Row_rev)) +
  #set up the platemap layout
      geom_point(data=expand.grid(seq(1, 12), seq(1, 8)), aes(x=Var1, y=Var2), color="grey90", fill="white", shape=21, size=8) +
  #Change the shape and colour of points for a variable
    geom_point(aes(colour = colour), size = 7) +
  #change the colours
    #scale_colour_viridis_c(guide = "colourbar", na.value="white") +
  # make into a plate plot
    theme_bdc_microtiter() +
  #fix the ratio of coordinates
    coord_fixed(ratio=(13/12)/(9/8), xlim=c(0.5, 12.5), ylim=c(0.6, 8.4)) +
  #add labels for the y axis
    scale_y_continuous(breaks=seq(1, 8), labels=LETTERS[8:1]) +
  #add labels for the x axis
    scale_x_continuous(position = "top", breaks=seq(1, 12)) +
  #Add a title
    labs(title=plate_names[i] , size = 6, colour = "QC Result") +
  #rotate legend guide because otherwise you can't see numbers:
    guides(fill = guide_colorbar(barwidth = 0.5, barheight = 10, title="Value"), colour = guide_legend(override.aes = list(size=4))) +
  #change the colours
   scale_colour_manual(values=c("red2", "green4")) +
  # make mimnimum point size bigger
    #scale_size_continuous(range = c(2,10)) +
  # make the font size in the legend bigger
    theme(legend.position="bottom", legend.text=element_text(size=10, colour = "black"), legend.title=element_text(size=10))
  
## print
plot_list[[i]] <- p
}

for (i in seq_along(plate_names)) {
    print(plot_list[[i]])
}

## load this library to use grid.arrange function
library(gridExtra)

## use this function to extract legend:
## source: https://stackoverflow.com/questions/13649473/add-a-common-legend-for-combined-ggplots
## source: https://github.com/hadley/ggplot2/wiki/Share-a-legend-between-two-ggplot2-graphs
g_legend<-function(a.gplot){
   tmp <- ggplot_gtable(ggplot_build(a.gplot))
   leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
   legend <- tmp$grobs[[leg]]
   return(legend)}

## get legend
mylegend<-g_legend(plot_list[[1]])

## make a final plot
QC_failed_plates <- grid.arrange(arrangeGrob(plot_list[[1]] + theme(legend.position="none") + labs(title = paste("Mutant 20", "\n", "plate 1")) + theme(plot.title = element_text(hjust = 0.5)),
                               plot_list[[2]] + theme(legend.position="none") + labs(title = paste("Mutant 20", "\n", "plate 2")) + theme(plot.title = element_text(hjust = 0.5)),
                               plot_list[[3]] + theme(legend.position="none") + labs(title = paste("Mutant 20", "\n", "plate 3")) + theme(plot.title = element_text(hjust = 0.5)),
                               plot_list[[4]] + theme(legend.position="none") + labs(title = paste("Mutant 20", "\n", "plate 4")) + theme(plot.title = element_text(hjust = 0.5)),
                               plot_list[[5]] + theme(legend.position="none") + labs(title = paste("Mutant 2", "\n", "plate 1")) + theme(plot.title = element_text(hjust = 0.5)),
                               plot_list[[6]] + theme(legend.position="none") + labs(title = paste("Mutant 2", "\n", "plate 2")) + theme(plot.title = element_text(hjust = 0.5)), nrow=2), 
                              mylegend, nrow=2,heights=c(10,2))

## plot
QC_failed_plates

TableGrob (2 x 1) "arrange": 2 grobs

## save
ggsave("/Users/Andy/GCSKO/GCSKO_analysis_git/images_to_export/QC_SS2_failed_plates.png", plot = QC_failed_plates, device = "png", path = NULL, scale = 1, width = 27, height = 25, units = "cm", dpi = 300, limitsize = TRUE)

Add bulk metadata

Add in bulk data predictions

hoo et al.

#Pb Prediction correlations with bulk data (asexual hoo): 

#Load in required package:
library(Hmisc)
#Cooerce expression data into a matrix and load in the reference timecourse data:
x10 <- as.matrix(ss2_mutants_final@assays$RNA@data)
rownames(x10) <- gsub("-", "_", rownames(x10))
#read in bulk data:
hoo<-as.matrix(read.table("/Users/Andy/GCSKO/GCSKO_analysis_git/data/Reference/hoo_berg2.txt",header=T, row.names=1))
#Make a blank dataframe in which to add prediction:
df <- data.frame(matrix(ncol = 4, nrow = 0))
colnames(df) <- c("Prediction(Spearman)","r(Spearman)","Prediction(Pearsons)","r(Pearsons)")
#Do correlations with bulk data using both Spearman and Pearson (and the top 1000 genes):
for (i in 1:ncol(x10))
{
  shared<-intersect(row.names(as.matrix(head(sort(x10[,i], decreasing=TRUE),1000))),row.names(hoo))
  step0<-rcorr(x10[shared,i],hoo[shared,1:12],type = "spearman")
  step1<-as.matrix(t(step0$r[2:13,1]))
  step2<-rcorr(x10[shared,i],hoo[shared,1:12],type = "pearson")
  step3<-as.matrix(t(step2$r[2:13,1]))
  step4<-cbind(colnames(step1)[which.max(step1)],step1[which.max(step1)],colnames(step3)[which.max(step3)],step3[which.max(step3)])
  colnames(step4) <- c("Prediction(Spearman)","r(Spearman)","Prediction(Pearsons)","r(Pearsons)")
  rownames(step4)<-colnames(x10)[i]
  df<-rbind(df,step4)
}
#Write out data into a csv file:
#write.csv(dfringr,file="/Users/ar19/Desktop/PhD/AR04_GCSKO_project/All_mutants_Feb_2018/predictionpbcombined.csv")
#Change the format of the output to make it more readable:
#gsub("Pb_","", dfringr[,1]) - Make predictions into 18hr.dat format:

#spearman:
df[,1] <- gsub("Pb_","", df[,1])
#Remove hr.dat from list:
df[,1] <- gsub("hr.dat","", df[,1])
#Check - dfringr[,1]
#Make into a number:
df[,1] <- as.numeric(df[,1])
df[,2] <- as.numeric(as.character(df[,2]))

#pearson:
df[,3] <- gsub("Pb_","", df[,3])
#Remove hr.dat from list:
df[,3] <- gsub("hr.dat","", df[,3])
#Check - dfringr[,1]
#Make into a number:
df[,3] <- as.numeric(df[,3])
df[,4] <- as.numeric(as.character(df[,4]))
#add to 10X object:
ss2_mutants_final <- AddMetaData(ss2_mutants_final, metadata = df)

Invalid name supplied, making object name syntactically valid. New object name is Prediction.Spearman.r.Spearman.Prediction.Pearsons.r.Pearsons.; see ?make.names for more details on syntax validity

Kasia’s data

Can also do with Kasia’s timecourse data:

kas<-as.matrix(read.table("/Users/Andy/GCSKO/GCSKO_analysis_git/data/Reference/AP2OETC.txt",header=T, row.names=1))
#Make a blank dataframe in which to add prediction:
dfs <- data.frame(matrix(ncol = 4, nrow = 0))
colnames(dfs) <- c("ID","Prediction","r (Pearson)")
#Do correlations with bulk data using both Spearman and Pearson (and the top 1000 genes):
for (i in 1:ncol(x10))
{
  shared<-intersect(row.names(as.matrix(head(sort(x10[,i], decreasing=TRUE),1000))),rownames(kas))
  step0<-rcorr(x10[shared,i],kas[shared,1:10],type = "spearman")
  step1<-as.matrix(t(step0$r[2:11,1]))
  step2<-rcorr(x10[shared,i],kas[shared,1:10],type = "pearson")
  step3<-as.matrix(t(step2$r[2:11,1]))
  step4<-cbind(colnames(step1)[which.max(step1)],step1[which.max(step1)],colnames(step3)[which.max(step3)],step3[which.max(step3)])
  colnames(step4) <- c("Prediction(Spearman)","r(Spearman)","Prediction(Pearsons)","r(Pearsons)")
  rownames(step4)<-colnames(x10)[i]
  dfs<-rbind(dfs,step4)
}
#Write out data into a csv file:
#write.csv(df,file="/Users/ar19/Desktop/PhD/AR04_GCSKO_project/All_mutants_Feb_2018/predictionkasiacombined.csv")

#Change the format of the output to make it more readable:
#gsub("Pb_","", dfs[,1]) - Make predictions into 18hr.dat format:
dfs[,1] <- gsub("X","", dfs[,1])
#Make into a number:
dfs[,1] <- as.numeric(dfs[,1])
#Make into a number:
dfs[,2] <- as.numeric(as.character(dfs[,2]))

#gsub("Pb_","", dfs[,1]) - Make predictions into 18hr.dat format:
dfs[,3] <- gsub("X","", dfs[,3])
#Make into a number:
dfs[,3] <- as.numeric(dfs[,3])
#dfs[,1]
#Make into a number:
dfs[,4] <- as.numeric(as.character(dfs[,4]))

colnames(dfs) <- c('Prediction(Spearman)_Kasia', 'r(Spearman)_Kasia', 'Prediction(Pearson)_Kasia', 'r(Pearson)_Kasia')
#add to Seurat:
#add to 10X object:
ss2_mutants_final <- AddMetaData(ss2_mutants_final, dfs)

Invalid name supplied, making object name syntactically valid. New object name is Prediction.Spearman._Kasiar.Spearman._KasiaPrediction.Pearson._Kasiar.Pearson._Kasia; see ?make.names for more details on syntax validity

Normalisation

normalise and find variable genes

## normalise
ss2_mutants_final <- NormalizeData(ss2_mutants_final, normalization.method = "LogNormalize", scale.factor = 10000)

Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

## find variable genes
ss2_mutants_final <- FindVariableFeatures(ss2_mutants_final, selection.method = "vst", nfeatures = 2000)

Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

## make a list of all genes in the dataset
all.genes <- rownames(ss2_mutants_final)

## scale data on all genes
ss2_mutants_final <- ScaleData(ss2_mutants_final, features = all.genes)

Centering and scaling data matrix

  |                                                                             
  |                                                                       |   0%
  |                                                                             
  |============                                                           |  17%
  |                                                                             
  |========================                                               |  33%
  |                                                                             
  |====================================                                   |  50%
  |                                                                             
  |===============================================                        |  67%
  |                                                                             
  |===========================================================            |  83%
  |                                                                             
  |=======================================================================| 100%

Dimensionality reduction

PCA

ss2_mutants_final <- RunPCA(ss2_mutants_final, features = VariableFeatures(object = ss2_mutants_final), verbose = FALSE)

PCA plot

DimPlot(ss2_mutants_final, reduction = "pca", group.by = "sub_identity_updated")

ElbowPlot(ss2_mutants_final, ndims = 30, reduction = "pca")

Find clusters

ss2_mutants_final <- FindNeighbors(ss2_mutants_final, dims = 1:21)

Computing nearest neighbor graph
Computing SNN

ss2_mutants_final <- FindClusters(ss2_mutants_final, resolution = 1)

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 2970
Number of edges: 97436

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.8634
Number of communities: 18
Elapsed time: 0 seconds

UMAP

Improve UMAP

#old way of running UMAP:
#ss2_mutants_final <-RunUMAP(ss2_mutants_final, reduction = "pca", dims = 1:10, n.neighbors = 150, seed.use = 1234, min.dist = 0.4, repulsion.strength = 0.03, local.connectivity = 150)
ss2_mutants_final <- RunUMAP(ss2_mutants_final, dims = 1:12)

20:09:17 UMAP embedding parameters a = 0.9922 b = 1.112
Found more than one class "dist" in cache; using the first, from namespace 'spam'
Also defined by ‘BiocGenerics’
20:09:17 Read 2970 rows and found 12 numeric columns
20:09:17 Using Annoy for neighbor search, n_neighbors = 30
Found more than one class "dist" in cache; using the first, from namespace 'spam'
Also defined by ‘BiocGenerics’
20:09:17 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:09:18 Writing NN index file to temp file /var/folders/wj/rztzclxn1t10cl2sk0plbf3r0000gn/T//RtmpXGhNV1/file43372c318018
20:09:18 Searching Annoy index using 1 thread, search_k = 3000
20:09:19 Annoy recall = 100%
20:09:21 Commencing smooth kNN distance calibration using 1 thread
20:09:23 Initializing from normalized Laplacian + noise
20:09:23 Commencing optimization for 500 epochs, with 117082 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:09:29 Optimization finished

p1 <- DimPlot(ss2_mutants_final, group.by = "ident", label = TRUE)
p2 <- DimPlot(ss2_mutants_final, group.by = "identity_updated")
p1 + p2

Cluster assignment

Interactive

plot <- FeaturePlot(ss2_mutants_final, features = "nFeature_RNA")   

HoverLocator(plot = plot, information = FetchData(ss2_mutants_final, vars = c("nFeature_RNA", "ident", "identity_updated")))

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be used`error_y.color` does not currently support multiple values.`error_x.color` does not currently support multiple values.`line.color` does not currently support multiple values.The titlefont attribute is deprecated. Use title = list(font = ...) instead.the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be used`error_y.color` does not currently support multiple values.`error_x.color` does not currently support multiple values.`line.color` does not currently support multiple values.The titlefont attribute is deprecated. Use title = list(font = ...) instead.

ngenes, nmito, ncounts visualisation

## add column to log counts so they are easier to visualise
ss2_mutants_final <- AddMetaData(ss2_mutants_final, log10(ss2_mutants_final@meta.data$nCount_RNA), col.name = "nCount_log10")

FeaturePlot(ss2_mutants_final, features = c("nCount_RNA", "nCount_log10", "nFeature_RNA", "percent.mt"))   

Violin plots cluster-by-cluster:

VlnPlot(object = ss2_mutants_final, features = "nFeature_RNA", pt.size = 0.01) +
  labs(x="Cluster",
       y="Genes per Cell",
       title = "Number of Genes per Cell") +
  theme(legend.position = "none")

VlnPlot(object = ss2_mutants_final, features = "nCount_log10", pt.size = 0.01) +
  labs(x="Cluster",
       y="Log10(Counts per Cell)",
       title = "Number of Counts per Cell") +
  theme(legend.position = "none")

VlnPlot(object = ss2_mutants_final, features = "percent.mt", pt.size = 0.01) +
  labs(x="Cluster",
       y="% Mitochondrial Reads",
       title = "Percentage of Mitochondrial Reads per Cell") +
  theme(legend.position = "none")

Expression of marker genes

# PBANKA-0515000 - p25 - female
# PBANKA-1319500 - CCP2 - female - used in 820 line
# PBANKA-1212600 - HAP2 - male
# PBANKA-0600600 - NEK3 - male
# PBANKA-1315700 - RON2 - (asexuals and some male?)
# "PBANKA-0416100" - dynenin heavy chain - male - used in 820 line
# PBANKA-1437500 - AP2-G - seuxal commitment gene
# PBANKA-0831000 - MSP1 - late asexual
# PBANKA-1102200 - MSP8 - early asexual (from Bozdech paper)
FeaturePlot(ss2_mutants_final, features = c("PBANKA-0515000", "PBANKA-1319500", "PBANKA-1212600","PBANKA-0600600", "PBANKA-1315700", "PBANKA-0416100", "PBANKA-1437500", "PBANKA-0831000", "PBANKA-1102200"), coord.fixed = TRUE)

Dotplot

This essentially allows us to view which genotypes sit in which clusters. Where we see that WT is not present, we can presume that this mutant has a unique transcriptome (either from the whole WT trajectory or it has arrested and so represents a WT developing form - we can find this out in the next stage of analaysis when we integrate with 10X data).

df_meta_data <- as.data.frame(ss2_mutants_final@meta.data)

dot_plot_df <- as.data.frame.matrix(table(df_meta_data$RNA_snn_res.1, df_meta_data$identity_updated))

dot_plot_df_pc <- (as.data.frame.matrix(prop.table(table(df_meta_data$RNA_snn_res.1, df_meta_data$identity_updated), margin = 2)) * 100)

dot_plot_df_pc$cluster <- rownames(dot_plot_df_pc)
library(dplyr)
dot_plot_df_pc_mutated <- mutate(dot_plot_df_pc)
library(reshape2)
dot_plot_df_pc_melted <- melt(dot_plot_df_pc, variable.name = "cluster")
colnames(dot_plot_df_pc_melted)[2] <- "identity"

df_meta_data <- as.data.frame(ss2_mutants_final@meta.data)

dot_plot_df <- as.data.frame.matrix(table(df_meta_data$RNA_snn_res.1, df_meta_data$identity_updated))

dot_plot_df_pc <- (as.data.frame.matrix(prop.table(table(df_meta_data$RNA_snn_res.1, df_meta_data$identity_updated), margin = 2)) * 100)

dot_plot_df_pc$cluster <- rownames(dot_plot_df_pc)
library(dplyr)
dot_plot_df_pc_mutated <- mutate(dot_plot_df_pc)
library(reshape2)
dot_plot_df_pc_melted <- melt(dot_plot_df_pc, variable.name = "cluster")

Using cluster as id variables

colnames(dot_plot_df_pc_melted)[2] <- "identity"

Mutants 3, 2, and 13 all seem to display this unique phenotype.

expression of 820 markers

library(ggplot2)
  p = ggplot(dot_plot_df_pc_melted,
             aes(y = factor(cluster),
                 x = factor(identity))) +
      geom_point(aes(colour=value, size=value)) + 
      scale_color_gradient(low="blue", high="red", limits=c( 1, max(dot_plot_df_pc_melted$value)), na.value="white") +
      theme_bw() +
      theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank())
  p = p +
      ylab("Cluster") +
      xlab("Identity") +
      theme(axis.text.x=element_text(size=12, face="italic", angle=45, hjust=1)) + 
      theme(axis.text.y=element_text(size=12, face="italic"))
  print(p)

Confirm life cycle designations using bulk correlations

## make plots 
plots <- FeaturePlot(ss2_mutants_final, features = c("PBANKA-1319500", "PBANKA-0416100"), blend = TRUE, combine = FALSE, coord.fixed = TRUE)

# Get just the co-expression plot, built-in legend is meaningless for this plot
#plots[[3]] + NoLegend()  

# Get just the key
#plots[[4]] 

# Stitch the co-expression and key plots together
plots[[3]] + NoLegend() + plots[[4]]/plot_spacer() + plot_layout(widths = c(2,1))

4. Cluster together with MCA data

Prepare

set up MCA data

FeaturePlot(ss2_mutants_final, features = c("Prediction.Spearman._Kasia", "Prediction.Spearman."))

set up MCA object

mca object

## read in data
counts_mca <- read.table("/Users/Andy/GCSKO/GCSKO_analysis_git/data/Reference/allcounts4.csv", header = TRUE, sep = ",", row.names=1, stringsAsFactors = TRUE)
dim(counts_mca)

[1] 5245 2432

anno_mca <- read.delim("/Users/Andy/GCSKO/GCSKO_analysis_git/data/Reference/allpheno8.2.csv", header = TRUE, sep = ",", row.names=1)
dim(anno_mca)

[1] 2432   49

## subset all blood stage cells
## find blood stages
blood_stages <- c("Merozoite", "Shz", "Male", "Schizont", "Female", "Trophozoite", "Ring")

blood_stage_cell_names <- rownames(anno_mca[anno_mca$ShortenedLifeStage2 %in% blood_stages, ])

## subset dataframes
counts_mca_blood_stage <- counts_mca[ ,colnames(counts_mca) %in% blood_stage_cell_names]
dim(counts_mca_blood_stage)

[1] 5245 1150

anno_mca_blood_stage <- anno_mca[rownames(anno_mca) %in% blood_stage_cell_names, ]
dim(anno_mca_blood_stage)

[1] 1150   49

## remove control cells:
non_control_cell_names <- rownames(anno_mca_blood_stage[anno_mca_blood_stage$Number_of_cells == 1, ])
counts_mca_blood_stage <- counts_mca_blood_stage[ ,colnames(counts_mca_blood_stage) %in% non_control_cell_names]
dim(counts_mca_blood_stage)

[1] 5245 1093

anno_mca_blood_stage <- anno_mca_blood_stage[rownames(anno_mca_blood_stage) %in% non_control_cell_names, ]
dim(anno_mca_blood_stage)

[1] 1093   49

## check
identical(rownames(counts_mca_blood_stage), rownames(counts))

[1] TRUE

n.b. a filter of 40 allows the merozoites to complete the ring in the asexuals

## make Seurat object with same filtering as main object because they are sequenced to same depth:
GCSKO_mca <- CreateSeuratObject(counts = counts_mca_blood_stage, meta.data = anno_mca_blood_stage, min.cells = 0, min.features = 120, project = "GCSKO")

Feature names cannot have underscores ('_'), replacing with dashes ('-')

GCSKO_mca

An object of class Seurat 
5245 features across 999 samples within 1 assay 
Active assay: RNA (5245 features, 0 variable features)

GCSKO_mca@meta.data$experiment <- "mca"
ss2_mutants_final@meta.data$experiment <- "ss2_mutants"

normalise object

VlnPlot(object = GCSKO_mca, features = "nFeature_RNA", pt.size = 0.01) +
  labs(x="Cluster",
       y="Genes per Cell",
       title = "Number of Genes per Cell") +
  theme(legend.position = "none") + geom_hline(yintercept=220)

VlnPlot(object = GCSKO_mca, features = "nCount_RNA", pt.size = 0.01) +
  labs(x="Cluster",
       y="Log10(Counts per Cell)",
       title = "Number of Counts per Cell") +
  theme(legend.position = "none")

Integrate

## normalise
GCSKO_mca <- NormalizeData(GCSKO_mca, normalization.method = "LogNormalize", scale.factor = 10000)

Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

## find variable genes
GCSKO_mca <- FindVariableFeatures(GCSKO_mca, selection.method = "vst", nfeatures = 2000)

Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

## make a list of all genes in the dataset
all.genes <- rownames(GCSKO_mca)
## scale data on all genes
GCSKO_mca <- ScaleData(GCSKO_mca, features = all.genes)

Centering and scaling data matrix

  |                                                                             
  |                                                                       |   0%
  |                                                                             
  |============                                                           |  17%
  |                                                                             
  |========================                                               |  33%
  |                                                                             
  |====================================                                   |  50%
  |                                                                             
  |===============================================                        |  67%
  |                                                                             
  |===========================================================            |  83%
  |                                                                             
  |=======================================================================| 100%

Analyse

Scale and PCA

## make list
ss2.mca.list <- list(GCSKO_mca, ss2_mutants_final)
ss2.mca.anchors <- FindIntegrationAnchors(object.list = ss2.mca.list, dims = 1:21, verbose = FALSE)

ss2.mca.integrated <- IntegrateData(anchorset = ss2.mca.anchors, dims = 1:21, verbose = FALSE, features.to.integrate = all.genes)

Adding a command log without an assay associated with it

Inspect PCS:

# set during IntegrateData
DefaultAssay(ss2.mca.integrated) <- "integrated"
# Run the standard workflow for visualization and clustering
ss2.mca.integrated <- ScaleData(ss2.mca.integrated, verbose = FALSE)
ss2.mca.integrated <- RunPCA(ss2.mca.integrated, npcs = 30, verbose = FALSE)

UMAP

## run UMAP
ss2.mca.integrated <- RunUMAP(ss2.mca.integrated, reduction = "pca", dims = 1:15)

20:11:28 UMAP embedding parameters a = 0.9922 b = 1.112
Found more than one class "dist" in cache; using the first, from namespace 'spam'
Also defined by ‘BiocGenerics’
20:11:28 Read 3969 rows and found 15 numeric columns
20:11:28 Using Annoy for neighbor search, n_neighbors = 30
Found more than one class "dist" in cache; using the first, from namespace 'spam'
Also defined by ‘BiocGenerics’
20:11:28 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:11:30 Writing NN index file to temp file /var/folders/wj/rztzclxn1t10cl2sk0plbf3r0000gn/T//RtmpXGhNV1/file4337275d142d
20:11:30 Searching Annoy index using 1 thread, search_k = 3000
20:11:31 Annoy recall = 100%
20:11:33 Commencing smooth kNN distance calibration using 1 thread
20:11:35 Initializing from normalized Laplacian + noise
20:11:36 Commencing optimization for 500 epochs, with 159194 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:11:43 Optimization finished

## add new column 
ss2.mca.integrated@meta.data$clusters_integrated <- ss2.mca.integrated@meta.data$seurat_clusters

ss2.mca.integrated@meta.data$clusters_integrated <- ifelse(is.na(ss2.mca.integrated@meta.data$clusters_integrated), ss2.mca.integrated@meta.data$ShortenedLifeStage2, ss2.mca.integrated@meta.data$clusters_integrated)

DimPlot(ss2.mca.integrated, pt.size = 0.01, label = TRUE, split.by = "experiment", group.by = "clusters_integrated") + coord_fixed() + theme_void()

ElbowPlot(ss2.mca.integrated, ndims = 30, reduction = "pca")

cluster

generate clusters:

VlnPlot(object = ss2.mca.integrated, features = "nFeature_RNA", pt.size = 0.01) + geom_hline(yintercept=100)

Save

save counts and pheno without anopheles and control cells

## copy old clusters
#ss2.mca.integrated <- AddMetaData(ss2.mca.integrated, ss2.mca.integrated@meta.data$RNA_snn_res.1, col.name = "pre_integration_clusters")

## generate new clusters
#ss2.mca.integrated <- FindNeighbors(ss2.mca.integrated, dims = 1:21)
#ss2.mca.integrated <- FindClusters(ss2.mca.integrated, resolution = 1, random.seed = 42, algorithm = 2)

save counts and pheno objects

write.csv(counts, file = "/Users/Andy/GCSKO/GCSKO_analysis_git/data_to_export/counts_2020_filtered.csv")
write.csv(pheno, file = "/Users/Andy/GCSKO/GCSKO_analysis_git/data_to_export/pheno_2020_filtered.csv")

save the session objects

pheno_filtered_cells <- ss2_mutants_final@meta.data
counts_filtered_cells <- ss2_mutants_final@assays$RNA@counts

write.csv(counts_filtered_cells, file = "/Users/Andy/GCSKO/GCSKO_analysis_git/data_to_export/counts_post_qc.csv")
write.csv(pheno_filtered_cells, file = "/Users/Andy/GCSKO/GCSKO_analysis_git/data_to_export/pheno_post_qc.csv")

Appendix

Session Info